Potential contribution of floral thermogenesis to cold adaptation, distribution pattern, and population structure of thermogenic and non/slightly thermogenic Symplocarpus species

Abstract The genus Symplocarpus in basal Araceae includes both thermogenic and non/slightly thermogenic species that prefer cold environments. If floral thermogenesis of Symplocarpus contributes to cold adaptation, it would be expected that thermogenic species have a larger habitat than non/slightly thermogenic species during an ice age, leading to increased genetic diversity in the current population. To address this question, potential distribution in past environment predicted by ecological niche modeling (ENM), genetic diversity, and population structure of chloroplast and genome‐wide single nucleotide polymorphisms were compared between thermogenic Symplocarpus renifolius and non/slightly thermogenic Symplocarpus nipponicus. ENM revealed that the distribution of S. nipponicus decreased, whereas that of S. renifolius expanded in the Last Glacial Maximum. Phylogeographic analyses have shown that the population structures of the two species were genetically segmented and that the genetic diversity of S. renifolius was higher than that of S. nipponicus. The phylogenetic relationship between chloroplast and nuclear DNA is topologically different in the two species, which may be due to the asymmetric gene flow ubiquitously observed in plants. The results of this study imply that floral thermogenesis of Symplocarpus contributes to expanding the distribution during an ice age, resulting in increased genetic diversity due to cold adaptation.


| INTRODUC TI ON
Thermogenic plants are defined as plants that can increase temperature in their reproductive organs (e.g., inflorescences, flowers, and cones) to at least 0.5°C above ambient temperature (Maekawa et al., 2022). Based on this definition, at least 90 plant species have been reported as thermogenic, half of which are gymnosperm cycads (Ito-Inaba et al., 2019;Tang, 1987), and the remaining half are angiosperms comprising mostly aroids (Araceae) and a few species of Nymphaeaceae, Magnoliacea, and Nelumbonaceae (Seymour, 2010).
Thermogenic aroids are independently observed in the subfamilies Orontioideae (basal Araceae with hermaphroditic flowers) and Aroideae (derived with monoecious flowers). Orontioideae includes only three genera (Symplocarpus, Lysichiton, and Orontium) that grow in northern temperate regions, such as eastern Asia and North America (Mayo et al., 1997). Among these, only Symplocarpus displays thermogenic features (Ito-Inaba et al., 2009). On the contrary, Aroideae is the major taxa of thermogenic aroids and contains many species with high heat-producing abilities. Floral thermogenesis is widely suggested to promote the dispersal of volatiles that may attract pollinators (Angioy et al., 2004;Gonçalves-Souza et al., 2017;Meeuse & Raskin, 1988), provide energetic benefits to adult insects (Seymour et al., 2003), facilitate fertilization (Li & Huang, 2009), and prevent freeze damage in plants (Knutson, 1974). However, evolutionary and ecological backgrounds remain unclear, such as how thermogenic plants have evolved adaptively or neutrally, why various plants have acquired thermogenesis in parallel, and if the benefits for plants are worth the energy consumption.
Among the thermogenic aroids, Symplocarpus is the only genus that favors inhabiting relatively cool temperate region. Several Symplocarpus species, such as American and Asian skunk cabbages (S. foetidus and S. renifolius, respectively), flowers upon thawing in early spring, maintain the spadix temperature at around 20°C despite below-freezing ambient air temperature for approximately 1 week (Knutson, 1974;Uemura et al., 1993). A previous study suggests the possibility that thermogenesis at low temperatures contributes to freezing avoidance by maintaining a high respiratory rate (Knutson, 1974). However, it remains largely unknown whether such floral thermogenesis contributes to survival or reproductive success in cold environment, particularly whether it is an adaptive trait. The markedly high heat-producing ability of Symplocarpus species is attributed to the increased cellular respiration facilitated by mitochondrial energy-dissipating proteins, such as alternative oxidase and uncoupling protein; however, this mechanism of heat generation comes at the cost of reduced ATP production (Ito-Inaba, Hida, Ichikawa, et al., 2008; Ito-Inaba, Hida, Mori, et al., 2008; Meeuse, 1975;Wagner et al., 2008). Therefore, given the high cost of floral thermogenesis, it can hardly be expected that the species will evolve entirely neutrally.
Currently, the genus Symplocarpus comprises six species (S. foetidus, S. renifolius, S. nipponicus, S. nabekuraensis, S. koreanus, and S. egorovii), half of which have been newly identified and described in the last 20 years (Lee et al., 2021;Otsuka et al., 2002;Pavlova & Nechaev, 2005). As mentioned above, S. foetidus and S. renifolius are well-known for their remarkable heat-producing abilities (Ito-Inaba et al., 2009;Knutson, 1974). Thermogenesis in the spadix has been previously reported in S. nabekuraensis and S. koreanus (Lee et al., 2021;Otsuka et al., 2011). Thermogenesis in S. nipponicus is rare or inexistent, and the individual size of this species is approximately one-third of S. renifolius (Otsuka et al., 2011). Meanwhile, thermogenesis in S. egorovii remains unknown. Among these species, S. renifolius, S. nabekuraensis, and S. nipponicus grow naturally in Japan and its surrounding areas. S. nipponicus and S. renifolius are widely distributed in these areas, whereas S. nabekuraensis is found only in higher altitude (Otsuka, 2002). From a phylogeographical perspective, it is of great interest that closely related species with different heat-producing abilities co-inhabit similar regions within Japan and its surrounding areas.
The common ancestor of the members of Simplocarpus diverged 8.8 million years ago and was independently exposed to past climate change (Lee et al., 2019). One clade underwent differentiation into thermogenic species (S. renifolius, S. nabekuraensis, and S. foetidus), and another clade underwent differentiation into non/slightly/ unknown-thermogenic species, S. nipponicus and S. egorovii. If floral thermogenesis contributes to cold adaptation in Symplocarpus, ecological and/or evolutionary consequences should be considered.
For example, if distributions between thermogenic and non/slightly thermogenic species in East Asia were different during an ice age, thermogenic species would expand to more northern and colder regions and distribute sporadically to suitable habitats or refugia.
Consequently, thermogenic species would have higher genetic divergence and clearly diverged population structures in the process where both species are distributed to the current areas.
To address this, we compared ecological niche modeling (ENM) predictions of potential distribution due to past climate change, niche divergence, genetic diversity, and population structure analysis of chloroplast DNA (cpDNA) and genome-wide single nucleotide polymorphisms (SNPs) between thermogenic S. renifolius and non/ slightly thermogenic S. nipponicus. S. nabekuraensis was omitted from the comparison since it is found in very limited areas and has a relatively small population. Our findings imply that floral thermogenesis contributes to the cold adaptation in Symplocarpus. and 2020 were dried using silica gel and stored at room temperature.

| Data and distribution modeling
Occurrence data for S. renifolius and S. nipponicus were obtained from the sampling locations and specimen information from the Tohoku University Botanical Gardens and downloaded from the Global Biodiversity Information Facility (GBIF.org, 12 November 2020) ( Figure 1). All records from museum specimens and GBIF were carefully checked against satellite images on Google Maps (https:// www.google.com/maps) and erroneous location data were removed.  Phillips & Dudík, 2008), Random Forest (RF;Breiman, 2001), and Support Vector Machine (SVM; Drake et al., 2006), to predict the potential current distribution of S. renifolius and S. nipponicus ( Figure S2). GLM and GAM are modeling methods based on regression analysis, while Maxent, RF, and SVM are based on machinelearning techniques. All models were fitted with five environmental variables before removing collinearities. In the Maxent parameter, the regularization multiplier beta was fitted from 1 to 10, with one as the best parameter, and the feature class was set by default automatic configuration. RF and SVM were tuned using tuneRF function of the randomForest package version 4.16-14 and tune.svm function of the e1071 package version 1.7-9 in R, respectively. The performance of these models was evaluated using the dismo package version 1.3-8 in R with five predictability indices: area under the receiver operating characteristic curve (AUC) (Swets, 1973), accuracy ([true positives + true negatives]/total n), sensitivity (true positives/ [true positives + false negatives]), specificity (true negatives/[false positives + true negatives]), and informedness (sensitivity + specificity − 1) ( Figure S2). The results showed that the best modeling methods were RF in S. renifolius and Maxent in S. nipponicus according to AUC, although none of the methods showed remarkably high performance. For interspecific comparison of the possible distribution range and factors determining the distribution, the Maxent method for each species was uniformly selected, but RF was not.
Maxent has been previously applied to closely related species, such as S. foetidus (Kim et al., 2018), and to species with small sample sizes and widespread distributions (Blair et al., 2013;Fuchs et al., 2018;van Proosdij et al., 2016). Assuming temporal stability of the ecological niche for both species, the constructed model was applied to LGM and MH climatic layers to predict the past distributions of both species.
Niche breadth and similarity between predicted niches were measured by Levin's niche breadth (Levins, 1968), Schoener's D (Schoener, 1968), and I (Warren et al., 2008) using ENMTools package, version 1.0.5 (Warren et al., 2021), in R. The differences in the niche between S. nipponicus and S. renifolius were evaluated by niche identity test and symmetric background test using ENMTools. The niche identity test addresses whether a pair of closely related species are effectively identical in their realized environmental distributions using null distribution constructed from actual occurrent points that are randomly reassigned. The symmetric background test addresses whether different environmental distributions reflect any underlying divergence in ecological tolerances or preferences using null distribution constructed from randomly chosen points from the broad region where both species live. For this test, random points around 500 km from the point of occurrence were used as background for two species. The null distributions were reconstructed from 999 random sampling replications and validated by env.D and env.I, which are improved indices of similarity between niches, with multiple corrections using false discovery rate. To investigate the association between bio19 and snowfall, the Mean Temperature of Coldest Quarter (bio11) of occurrence locations was obtained from the current environments using raster package (version 3.5-21) and tested using the Wilcoxon rank sum test with the coin package (version 1.4-2) in R.

| DNA extraction and sequencing
Total genomic DNA was extracted using hexadecyltrimethylammonium bromide (CTAB) method (Murray & Thompson, 1980), with minor modification. To acquire haplotypes of cpDNA, primers for rbcL (Kress et al., 2009) and psbA-trnH intergenic spacer (Sang et al., 1997) were used for 1st PCR. These primers were preceded by one to three random bases (N, NN, or NNN). For the second PCR, complementary sequences for the binding sites of the Illumina sequencing flow cell and barcode sequences for each sample were added to the first PCR products. Both PCRs were performed using PrimeSTAR HS DNA Polymerase (Takara). These templates were pooled and sequenced using an Illumina MiSeq sequencer (Illumina) with Reagent Kit v2 (500 cycle, paired-end).
To acquire genome-wide SNPs as nuclear markers, multiplexed inter-simple sequence repeats genotyping by sequencing (MIGseq) was performed (Suyama et al., 2022). MIG-seq is a reducedrepresentation DNA sequencing method that uses a high-throughput sequencing platform to detect putatively neutral genome-wide SNPs adjacent to microsatellite regions. It is used for low-quality DNA, such as museum specimens, because of the PCR amplification-based approach (Strijk et al., 2020). The MIG-seq library was sequenced using an Illumina MiSeq sequencer with a Reagent Kit v3 (150 cycle, paired-end).

| Genotyping
Reads obtained from cpDNA sequencing were analyzed according to Claident protocol (Tanabe & Toju, 2013) using the recommended parameters (https://www.claid ent.org). Raw data obtained from the sequencer were called in fastq format using bcl2fastq v2 (Illumina), F I G U R E 1 Occurrence sites of Symplocarpus renifolius and Symplocarpus nipponicus in Japan, Sakhalin, and Northeast Asia obtained from GBIF database, museum specimen, and sampling points for this study. Black dots indicate background sites using ecological niche modeling. and the reads were demultiplexed. Because of different insert sizes, paired-end reads of psbA-trnH were merged at overlapping regions, and rbcL reads were concatenated. Low-quality or noisy sequence were filtered out, and passing sequences were clustered by 100% identity. The most frequently acquired reads of each sample were obtained as representative sequences. The representative sequences of all samples accounted for more than 50% of the obtained sequences, and the genera were confirmed using the Basic Local Alignment Search Tool (BLAST; https://blast.ncbi.nlm.nih.gov/ Blast.cgi). These treatments were applied to reduce contaminations.
The representative sequences were aligned using MAFFT (Katoh & Standley, 2013), and 5-bp around gaps were trimmed for proper alignment. Lastly, the representative sequences of each sample were concatenated.
MIG-seq reads were trimmed by removing adapter sequences and low-quality reads using Trimmomatic version 0.39 (Bolger et al., 2014). The filtered reads were genotyped using stacks version 2.4 (Rochette et al., 2019) for de novo assembly. The parameter of minimum depth of coverage required to classify reads as "stack" was set to three (m = 3). The number of mismatches allowed between stacks and to align secondary reads were set to one and two, respectively (M = 1, n = 1, N = 2). These mismatch parameters are relatively low to avoid misassembly of paralogs since S. renifolius is tetraploid. This reduces the number of SNPs but is more conservative. The stacks were genotyped as SNPs if they were typed as at least 70% of the target species, minor alleles were counted to more than two, and observed heterozygosity was less than 60% using the "populations" command. The parameter settings were increased to obtain more SNPs (M = 2, n = 2, N = 4) yielded qualitatively similar results.
Populations with four or more amplified individuals were used for nuclear DNA analysis. A total of 135 individuals from 15 populations of S. renifolius (mean of 9.0 individuals per population, 5-24) and 92 individuals from eight populations of S. nipponicus (mean of 11.5 individuals per population, 6-24) were analyzed ( Table S1). The rate of the polymorphic loci, nucleotide diversity (π), observed heterozygosity (H o ), and the number of private alleles were calculated by the "populations" command in Stacks as summary statistics of genetic diversity for each population of both species. Phylogenetic relationships among individuals were determined by maximum likelihood estimation using RAxML version 8.2.12 (Stamatakis, 2014) with GTRGAMMA model and visualized using ggtree (Yu et al., 2017).
The population structure of each species was determined using ADMIXTURE version 1.3.0 (Alexander & Lange, 2011). The numbers of clusters were tested from 1 to 15, and the other parameters were set by default values. The best K values were selected based on cross-validation error.

| Comparison of potential distribution and response to past climate change
To investigate the relationship between plant thermogenesis and cold adaptation, the potential distributions of thermogenic S. renifolius and non/slightly thermogenic S. nipponicus were determined and compared using Maxent. Maxent models for both species showed high accuracy in predicting the observed distributions, with high AUC values (>0.9) for the best parameters LGM (Figure 3a). The similarity of predicted distribution between two species in the LGM was lower than that in the current and MH in both indexes (Figure 3b). Comparing the similarity between the ages, the similarity between the MH and the LGM in S. nipponicus was lower than that in S. renifolius in both indexes (Figure 3c), suggesting different distribution shift between species from the LGM to the current environment. These notions were further confirmed by the visualization (Figure 3d,e). The distribution of S. nipponicus in the MH was similar to the current distribution, while the distribution with higher probability in the LGM was narrower and more sporadic (Figure 3d). In contrast, when compared to the current distribution, the distribution of S. renifolius in the MH was narrower and more sporadic, with a lower probability, while that in the LGM was broader (Figure 3e; Figure S4). Based on the Maxent model, the potential distribution of S. renifolius was widely expanded to the colder and snowier northern region than that of S. nipponicus in the Japanese islands and Sakhalin, which roughly corresponds to its current occurrence (Figures 1 and 3). Overall, these results showed that the niche of S. nipponicus has kept the breadth and migrated since the LGM, the niche of S. renifolius has decreased since the LGM, and the snowfall in the coldest quarter contributed the most to the distribution pattern of the two species.

| Phylogenetic and geographic relationship of genome-wide SNPs by MIG-seq
To perform a detailed comparison of phylogenetic relationships and population structure within the populations of S. nipponicus and S. renifolius by nuclear genome sequencing, MIG-seq was performed using a high-throughput sequencer ( Figure 5). MIG-seq is a powerful tool for genome-wide identification of SNPs by sequencing the region between microsatellites. The evolutionary timescale in nuclear genome is largely different from that in the chloroplast genome, and generally, evolutionary processes in the nuclear genome are more rapid than those in the chloroplast genome (Drouin et al., 2008) (Figure 5a; Figure S5). The percentage of polymorphic loci, nucleotide diversity, and observed heterozygosity in the nuclear genome of S. renifolius was significantly higher than those in S. nipponicus (p < .005, Welch's t-test with Benjamini and Hochberg correction (Benjamini & Hochberg, 1995) for multiple testing), but the number of the private allele was not significant (p = .208) (Figure 5b; Table S1). In S. renifolius with higher genetic diversity, the phylogenetic relationship was topologically inconsistent with the haplotype network analyses using cpDNA (Figures 4a and 5c).

| Contribution of floral thermogenesis to niche divergence, cold adaptation, and pollination ecology
This study shows that the species exhibiting higher thermogenesis also lives in colder climates. The current distributions between species were significantly different from those predicted from random using ENM comparison, suggesting niche divergence. The different responses of the distributions between species were illustrated by the past fluctuated environment of the MH, which was hotter than the current environment, and the LGM, which  (Table S1). Asterisk indicates a significant difference between species using Welche's t-test with Benjamini and Hochberg correction. (c) Phylogenetic tree using RAxML. Dot colors indicate haplotype of cpDNA ( Figure 4) and colored bars around the trees indicate population structure of best K in (a). If floral thermogenesis contributes to cold adaptation in the thermogenic Symplocarpus, genetic vestiges remain, such that the genetic diversity of thermogenic species is higher than that of non/ slightly thermogenic species because of the avoidance of distribution reduction in an ice age. Comparative analysis between thermogenic S. renifolius and non/slightly thermogenic S. nipponicus has revealed the different genetic divergence in cpDNA and nuclear DNA, as well as the topological discordance of phylogenetic relationship between chloroplast and nuclear genome in S. renifolius and fragmentary population genetic structure in nuclear DNA (Figures 4 and 5).
Topological discordance between organelle and nuclear genome has occasionally been observed in animals (Ford et al., 2019;Singhal & Moritz, 2012) and plants (Hirota et al., 2021;Petit et al., 2005;Xu et al., 2021). This may be due to the different patterns of gene flow between maternally inherited organelles and biparentally inherited nuclei, and sex-biased dispersal between seed and pollen. Although seeds usually fall into muddy substrates or are carried off by animals and floods, almost all seeds of both species are consumed by rodents in Japan (Otsuka & Kitano, 2003) and the range of seed dispersal by rodents has been estimated to be approximately 10 m in S. renifolius (Abe et al., 2006;Wada & Uemura, 1994). In contrast, pollen are that thermogenic species will not be able to attract pollinators and provide benefits to them in early spring. This suggestion applies to thermogenic species that prefer cold environments, although climatic change will also affect to other thermogenic plants according to role of heat in each species.
Moreover, a comparison between thermogenic and non/slightly thermogenic species revealed different response to the amount of snowfall and different historical distributions using Maxent. It also revealed different ranges of pollen dispersal from phylogeographic discordance between cpDNA and genome-wide SNPs. Our findings suggest that floral thermogenesis is a consequence of adaptive evolution due to cold environment and pollination, but it remains a possibility that the results are caused by the history of distribution formation after speciation, instead of adaptation. The possibility that phylogeographic discordance is caused by positive selection could not be investigated in this study because of limited genetic information. Closely related species available for verification are L. camtschatcensis, which is non-thermogenic and broadly distributed in Japan and Russia, and S. nabekuraensis, which is thermogenic and locally distributed in high-altitude areas in Japan, but an endangered species. To elucidate the adaptive signature of thermogenesis in Symplocarpus, genome sequencing, understanding thermogenic mechanisms, and detecting the natural selection of thermogenic genes are necessary.

| Phylogeography and speciation process in Symplocarpus
The phylogenetic relationship between cpDNA and the nuclear genome of Symplocarpus was consistent with previous studies, except for the outgroup L. camtschatcensis (Kitano et al., 2005;Lee et al., 2019;Nie et al., 2006).  (Lee et al., 2021;Pavlova & Nechaev, 2005). According to the habitat ranges estimated by occurrence data in Japan using Maxent model, the potential current distribution of S. renifolius in